#coding:utf8
'''
Created on 2018年8月27日

@author: Administrator
'''
from pylab import *

# Set up initial condition
x = -10.0
y = -7.0
z = 35.0
# Horizonal and vertical 2D plot projections
xt = []
yt = []

dt = 0.001
sigma = 10.0
R = 28.0
b = -8.0/3.0
theta = 3.0 * pi / 4.0  # Rotate about z-axis

for t in arange(0,10,dt):
    dxdt = sigma * (-x + y)
    dydt = R*x - x*z - y
    dzdt = b*z + x*y
    x = x + dxdt*dt
    y = y + dydt*dt
    z = z + dzdt*dt
    # Rotate orientation in 3D for a better view
    #   about z-axis
    x_rot = cos(theta) * x - sin(theta) * y
    y_rot = sin(theta) * x + cos(theta) * y
    xt.append(x_rot)
    yt.append(z)

# Setup the plot
xlabel('x(t)+y(t)') # set x-axis label
ylabel('z(t)') # set y-axis label
# Set plot title
PlotTitle = 'Lorenz ODE with (sigma,b,R) = (%5.2f,%5.2f,%5.2f)' % (sigma,b,R)
title(PlotTitle)
# Plot the time series
plot(xt,yt, 'k-')
# Make sure the orbit appears in a unit square
axis('equal')
axis([-30.0, 30.0, 0.0, 60.0])

# Use command below to save figure
#savefig('LorenzOrbit', dpi=600)

# Display the plot in a window
# show()

print xt[0:50]
print yt[0:50]


